function [params_hat, covMat, params_ses, modelFits] = estimate_model_by_blocks(modelId, data, NumBlocks, varargin)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function estimates the model by maximum likelihood (using Newton-Raphson algorithm).
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% modelId;					integer
	% data:			           object:
	% 	.dims:				     	object:
	%		.dims1:						cell(NumParts,1)
	%			{ii}:						integer: gives dim1_ii that corresponds to Xparts{ii}
	%		.Xpart_2_NumX:				1 x NumParts: gives integers (dim2 of Xparts{ii}.X)
	%		.Xpart_2_NumX_FEs:			1 x NumParts: gives integers (dim2 of Xparts{ii}.X_FEs)
	%		.Xpart_2_Num_FE_vals:		cell(1,NumParts)
	%			{ii}:						1 x NumX_FEs_i  --> gives integer: number of possible values for Xparts{ii}.X_FEs(:,ff)
	%		.NumFEvals2Keep:			cell(1,NumParts) --> gives integer
	%		.NumParts
	%		.NumObs
	%		.NumParams
	%    .LL_offset:               1 x 1
	%	 .Xparts:				   cell(NumParts,1)
	%	 	{ii}:				       object
	%			.X:					       dim1_ii x NumX
	%			.X_FEs:					   dim1_ii x NumX_FEs
	%			.NumX_FE_vals:			   1 x NumX_FEs: gives integer
	%    .N:                       NumObs x 1    (if modelId==1)
	%    .logLambdaOffset:         NumObs x 1     (if modelId==2)
	%    .Y:                       NumObs x 1
	% varargin{1}=params0:		   NumParams x 1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% params_hat:				   NumParams x 1
	% covMat:					   NumParams x NumParams
	% params_ses:				   NumParams x 1
	% modelFits:				   object:
	%	.LL:							1 x 1
	%	.AIC:							1 x 1
	%	.BIC:							1 x 1
	%	.MAE:							1 x 1
	%	.MSE:							1 x 1
	%	.RMSE:							1 x 1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	% Read dimensions
	NumParams = data.dims.NumParams;
	NumObs    = size(data.Y,1);
	disp(sprintf('Number of observations: %d', NumObs));
	disp(sprintf('Number of parameters: %d',   NumParams));
	
	% Read params or set default value
	if nargin >= 4
		params0 = varargin{1};
	else
		params0 = zeros(NumParams,1);
	end
		
	%%% Define objective function
	if data.dims.dimType ~= 2 && data.dims.dimType ~= 3; error('This below assumes dimType=2 or dimType=3'); end;
	% Read LL_offset
	LL_offset = data.LL_offset;
	
	% Read dimensions
	dims = data.dims;
	NumParts = dims.NumParts;
	if dims.dimType == 2
		T = dims.dims1{2};
		K = dims.dims1{3};
	else if dims.dimType == 3
		T = dims.dims1{4};
		K1 = dims.dims1{5};
		K2 = dims.dims1{6};
	end; end;
	
	%%% Split the data by blocks
	NumTs_perBlock = ceil(T/NumBlocks);
	tt_offset = 0;
	data_blocks = cell(NumBlocks, 1);
	
	for bb = 1:NumBlocks
		block_tts = tt_offset+1:min(tt_offset+NumTs_perBlock, T);
		T_bb = length(block_tts);
	
		%%% Make Xparts_bb and dims_bb
		%% Pieces common across all bb
		dims_bb.dimType             = dims.dimType;
		dims_bb.Xpart_2_NumX        = dims.Xpart_2_NumX;
		dims_bb.Xpart_2_NumX_FEs    = dims.Xpart_2_NumX_FEs;
		dims_bb.Xpart_2_Num_FE_vals = dims.Xpart_2_Num_FE_vals;
		dims_bb.NumFEvals2Keep      = dims.NumFEvals2Keep;
		dims_bb.NumParts            = dims.NumParts;
		dims_bb.NumParams           = dims.NumParams;
		
		if dims.dimType == 2
			Xparts_bb = cell(3,1);
			Xparts_bb{3} = data.Xparts{3}; % K x 1
			
			%% Pieces that depend on bb
			% Xpart 1
			NumX_part1     = size(data.Xparts{1}.X, 2);
			NumX_FEs_part1 = size(data.Xparts{1}.X_FEs, 2);
			Xpart1_X     = reshape(data.Xparts{1}.X,     [T K NumX_part1]);
			Xpart1_X_FEs = reshape(data.Xparts{1}.X_FEs, [T K NumX_FEs_part1]);
			Xpart1_X     = reshape(Xpart1_X(block_tts,:,:),     [T_bb*K NumX_part1]);     % (T_bb*K) x NumX
			Xpart1_X_FEs = reshape(Xpart1_X_FEs(block_tts,:,:), [T_bb*K NumX_FEs_part1]); % (T_bb*K) x NumX_FEs_part1
			Xpart1.X            = Xpart1_X;
			Xpart1.X_FEs        = Xpart1_X_FEs;
			Xpart1.NumX_FE_vals = data.Xparts{1}.NumX_FE_vals;
			Xparts_bb{1} = Xpart1; 
			clear Xpart1 NumX_part1 NumX_FEs_part1 Xpart1_X Xpart1_X_FEs;
			
			% Xpart 2
			Xpart2.X            = data.Xparts{2}.X(block_tts,:); % T_bb x NumX
			Xpart2.X_FEs        = data.Xparts{2}.X_FEs(block_tts,:); % T_bb x NumX_FEs
			Xpart2.NumX_FE_vals = data.Xparts{2}.NumX_FE_vals;
			Xparts_bb{2} = Xpart2; 
			clear Xpart2;
			
			% dims_bb.NumObs
			dims_bb.NumObs = T_bb*K;
			
			% Make Y, N and/or logLambdaOffset
			Y_bb = reshape(data.Y, [T K]); % T x K
			Y_bb = reshape(Y_bb(block_tts,:), [T_bb*K 1]); % (T_bb*K) x 1
			data_bb.Y = Y_bb;
			if isfield(data, 'N');
				N_bb = reshape(data.N, [T K]); % T x K
				N_bb = reshape(N_bb(block_tts,:,:), [T_bb*K 1]); % (T_bb*K) x 1
				data_bb.N = N_bb;
			end;
			if isfield(data, 'logLambdaOffset');
				logLambdaOffset_bb = reshape(data.logLambdaOffset, [T K]); % T x K
				data_bb.logLambdaOffset = logLambdaOffset_bb(block_tts,:,:); % T_bb x K
			end;
			
		else if dims.dimType == 3
		
			Xparts_bb = cell(6,1);
			Xparts_bb{1} = data.Xparts{1}; % (K1*K2) x 1
			Xparts_bb{5} = data.Xparts{5}; % K1 x 1
			Xparts_bb{6} = data.Xparts{6}; % K2 x 1
			
			%% Pieces that depend on bb
			% Xpart 2
			NumX_part2     = size(data.Xparts{2}.X, 2);
			NumX_FEs_part2 = size(data.Xparts{2}.X_FEs, 2);
			Xpart2_X     = reshape(data.Xparts{2}.X,     [T K1 NumX_part2]);
			Xpart2_X_FEs = reshape(data.Xparts{2}.X_FEs, [T K1 NumX_FEs_part2]);
			Xpart2_X     = reshape(Xpart2_X(block_tts,:,:),     [T_bb*K1 NumX_part2]);     % (T_bb*K1) x NumX
			Xpart2_X_FEs = reshape(Xpart2_X_FEs(block_tts,:,:), [T_bb*K1 NumX_FEs_part2]); % (T_bb*K1) x NumX_FEs_part2
			Xpart2.X            = Xpart2_X;
			Xpart2.X_FEs        = Xpart2_X_FEs;
			Xpart2.NumX_FE_vals = data.Xparts{2}.NumX_FE_vals;
			Xparts_bb{2} = Xpart2; 
			clear Xpart2 NumX_part2 NumX_FEs_part2 Xpart2_X Xpart2_X_FEs;
			
			% Xpart 3
			NumX_part3     = size(data.Xparts{3}.X, 2);
			NumX_FEs_part3 = size(data.Xparts{3}.X_FEs, 2);
			Xpart3_X     = reshape(data.Xparts{3}.X,     [T K2 NumX_part3]);
			Xpart3_X_FEs = reshape(data.Xparts{3}.X_FEs, [T K2 NumX_FEs_part3]);
			Xpart3_X     = reshape(Xpart3_X(block_tts,:,:),     [T_bb*K2 NumX_part3]);     % (T_bb*K2) x NumX
			Xpart3_X_FEs = reshape(Xpart3_X_FEs(block_tts,:,:), [T_bb*K2 NumX_FEs_part3]); % (T_bb*K2) x NumX_FEs_part3
			Xpart3.X            = Xpart3_X;
			Xpart3.X_FEs        = Xpart3_X_FEs;
			Xpart3.NumX_FE_vals = data.Xparts{3}.NumX_FE_vals;
			Xparts_bb{3} = Xpart3; 
			clear Xpart3 NumX_part3 NumX_FEs_part3 Xpart3_X Xpart3_X_FEs;
			
			% Xpart 4
			Xpart4.X            = data.Xparts{4}.X(block_tts,:); % T_bb x NumX
			Xpart4.X_FEs        = data.Xparts{4}.X_FEs(block_tts,:); % T_bb x NumX_FEs
			Xpart4.NumX_FE_vals = data.Xparts{4}.NumX_FE_vals;
			Xparts_bb{4} = Xpart4; 
			clear Xpart4;
			
			% dims_bb.NumObs
			dims_bb.NumObs              = T_bb*K1*K2;
			
			% Make Y, N and/or logLambdaOffset
			Y_bb = reshape(data.Y, [T K1*K2]); % T x (K1*K2)
			Y_bb = reshape(Y_bb(block_tts,:), [T_bb*K1*K2 1]); % (T_bb*K1*K2) x 1
			data_bb.Y = Y_bb;
			if isfield(data, 'N');
				N_bb = reshape(data.N, [T K1 K2]); % T x K1 x K2
				N_bb = reshape(N_bb(block_tts,:,:), [T_bb*K1*K2 1]); % (T_bb*K1*K2) x 1
				data_bb.N = N_bb;
			end;
			if isfield(data, 'logLambdaOffset_dim1');
				data_bb.logLambdaOffset_dim1     = data.logLambdaOffset_dim1(block_tts,:); % T_bb x K1
				data_bb.logLambdaOffset_dim2 = data.logLambdaOffset_dim2(block_tts,:,:); % T_bb x 1 x K2
			end;
		end; end;
		
		if isfield(data, 'spellDurations');
			data_bb.spellDurations = data.spellDurations(block_tts,:); % T_bb x 1
		end;
		
		for pp = 1:NumParts
			dims_bb.dims1{pp}    = size(Xparts_bb{pp}.X,1);
		end
		
		% Save everything
		data_bb.dims   = dims_bb;
		data_bb.Xparts = Xparts_bb;
		data_bb.LL_offset = 0;
		data_blocks{bb} = data_bb;
		
		% Update tt_offset
		tt_offset = tt_offset + T_bb;
	end
	clear data;
	data.LL_offset = LL_offset;
	data.blocks = data_blocks;
	
	myfunc = @(params_val) loglikelihood_by_blocks(modelId, data, params_val, true);
	
	%%% Estimate by Newton-Raphson
	tic
	params_hat = Newton_Raphson(myfunc, params0, 1e-8, 1000, 1, true);
	toc
	
	%%% Calculate standard errors
	[LL,~,FisherInfo] = myfunc(params_hat);
	LL = -LL;
	covMat = inv(FisherInfo);
	params_ses = sqrt(diag(covMat)); % Standard errors
	
	%%% Compute AIC and BIC
	modelFits.LL = LL;
	modelFits.AIC = 2*NumParams - 2*LL;
	modelFits.BIC = NumParams*log(NumObs) - 2*LL;
	
	modelFits.MAE = -1;
	modelFits.MSE = -1;
	modelFits.RMSE = -1;

end
